Linear instability of planar shear banded flow 
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We study the linear stability of planar shear banded flow with respect to perturbations with 
wavevector in the plane of the banding interface, within the non local Johnson Segalman model. We 
find that perturbations grow in time, over a range of wavevectors, rendering the interface linearly 
unstable. Results for the unstable eigenfunction are used to discuss the nature of the instability. 
We also comment on the stability of phase separated domains to shear flow in model H. 

PACS numbers: 47.50. +d Non-Newtonian fluid flows- 83.60.Wc Flow instabilities - 47.55. Kf Multi-phase 
and particle laden flows - 61.25.Hq Macromolecular and polymer solutions; polymer melts; swelling 



Complex fluids such as wormlike micellar surfac- 
tants l], lamellar onion phases @, polymer solutions Q 
and soft glasses commonly undergo flow instabilities 
and flow-induced transitions that result in spatially het- 
erogeneous "shear banded" states. This effect is captured 
by several notable rheological models [|| in which the un- 
derlying constitutive curve of shear stress vs. shear rate, 
1^(7), is non-monotonic (Fig. allowing the coexis- 
tence of bands of differing shear rate at common shear 
stress, Fig. [21 However, most theoretical studies have 
considered only one spatial dimension (ID) 0,01 , normal 
to the interface between the bands (the flow gradient di- 
rection, y). The stability of ID banded profiles in higher 
dimensions has been implicitly assumed, but is in fact 
an open question. In this Letter, therefore, we study nu- 
merically the linear stability of ID planar shear banded 
profiles with respect to perturbations with wavevectors 
in the interfacial plane (x, z) — (flow, vorticity). 

We work within the Johnson Segalman (JS) model 
modified to include non local diffusive terms 0. These 
account for gradients in the order parameters across the 
banding interface, conferring a surface tension. This 
"dJS" model is often taken as a paradigm of shear band- 
ing systems. Our main result will be that interfacial fluc- 
tuations typically grow in time, rendering the ID banded 
profile linearly unstable. This potentially opens the way 
to non trivial interfacial dynamics and could form a start- 
ing point for understanding an emerging body of data 
revealing erratic fluctuations of shear banded flows [Io| . 
This work is a timely counterpart to new techniques for 
measuring interfacial dynamics [Tlj . It is also relevant 
industrially, to processing instability and oil extraction. 

The model is defined as follows. The generalised 
Navier Stokes equation for a viscoelastic material in a 
Newtonian solvent of viscosity 77 and density p is: 

/0 (d t + V.V)V = V-CE + ryVV-PI), (1) 

where V(R) is the velocity field and £(R) the viscoelas- 
tic part of the stress. For homogeneous planar shear, 
V = 2/7X, the total shear stress T xy — £,^(7) + 777. The 
pressure P is determined by incompressibility, 

V.V = 0. (2) 



The viscoelastic stress evolves with dJS dynamics ja, || 
X I 2 

S= 2GD - - + -V 2 S, (3) 

T T 

with plateau modulus G and relaxation time r. The non 
local diffusive term accounts for spatial gradients across 
the interface between the bands. It arises naturally in 
models of liquid crystals, and diffusion of strained poly- 
mer molecules 12]. The time derivative 

£= (d t + V • V) £ - a(D ■ S + £ • D) - (S • ft - SI - £), 

in which D and are the symmetric and antisymmetric 
parts of the velocity gradient tensor, (VV) Q ^ = d a vp. 
The "slip parameter" a measures the non-affinity of de- 
formation of the viscoelastic component |8|]. Slip occurs 
for I a I < 1. The underlying constitutive curve T xy {^) is 
then capable of the non-monotonic behaviour of Fig. 

Within this model we consider planar shear between 
infinite, flat parallel plates at y = 0, L. We use units in 
which G = 1, r = 1 and L = 1; and boundary conditions 
at y = 0, 1 of d y T, a p — V a, f3 for the viscoelastic stress, 
with no slip and no penetration for the velocity. 

For an imposed shear rate 7 in the region of decreasing 
stress, dTxy/d^ < 0, homogeneous flow is unstable [l3| . 
A ID analysis in the flow gradient dimension then pre- 
dicts a separation into two bands of differing shear rates 
71,72 at common shear stress, T^, separated by an in- 
terface of width 0(1). As the applied shear rate 7 is 
tracked across the banding regime, the relative width- 
fraction of the bands adjusts to maintain the constraint 
J dy A f(y) — 7, while 71,72 and lb stay constant, leading 
to a plateau in the steady state flow curve (Fig.0. 

We verified this ID scenario by numerically evolving 
Eqns.^to|31 allowing spatial variations only in the flow- 
gradient direction y. We used a Crank Nicholson algo- 
rithm |bl| within a finite difference scheme on a uniform 
mesh of "full" points 7/0, yi---yN b ^ al . for £ and staggered 
"half" points yi , ys ...y Nb ^_ 1 for V. We evolved with 
time-step Dt for a time /j max to steady state, checking for 
convergence to the limit A^asc — * 00, Dt — > 0, /j max — > 00. 

The resulting flow curve is shown in Fig. ^ A typical 
steady state shear banded profile V(y), £(?/) is given 
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FIG. 1: Underlying constitutive curve; steady state flow 
curve, a = 0.3, r/ = 0.05. Banding occurs on the plateau. 
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FIG. 2: ID banded profile, with spatial gradients restricted 
to the flow gradient direction, y. 7 = 2.0, towards the left of 
the plateau in Fig. [I] I = 0.01, N hase = 800. 



in Fig. [21 The velocity normal to the interface V y = 
in this ID profile. The smooth variation of the order 
parameters across the interface results from the diffusive 
term in Eqn. which confers an interface width 0(1). 
This is in contrast to local models (I = 0) in which the 
interface is a sharp discontinuity. In fact, local models 
are pathological in the sense that the banded state is not 
uniquely selected, but depends on flow history 

The linear stability of the sharply banded profiles 
of local models was studied by previous authors. Re- 
nardy |15| found instability with respect to interfacial 
fluctuations of high wavevector, q x — > 00, in the local JS 
model restricted to the case of a thin high shear band. 
McLeish [l6j studied capillary flow, for general band 
thickness. He demonstrated a long wavelength (q x — » 0) 
instability due to the jump in normal stresses across the 
interface. This mechanism was also discussed in Ref. 

Here we study numerically the non local case, in which 
the ID banded profile is uniquely selected 0, 0- The 
non zero interfacial width, I <C L, now confers a surface 
tension, which was absent from the local case. We study 
general band thicknesses and the full (velocity, vorticity) 
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FIG. 3: Real part of the eigenvalue of the most unstable mode. 
a = 0.3, n = 0.05, 7 = 2.0, Reynolds number p/rj = 0. The 
data for I — 0.01 correspond to the base profile in Fig. [5] 
Symbols: data. Solid lines: cubic splines. 



plane of perturbation wavevectors (q Xl q z ). 

We linearised the model equations to for small per- 
turbations (lower case) about the (upper case) base pro- 
file, i>(x, y, z, t) = + 4>q(y) exp(w q t + iq x x + iq z z). 
The vector $ comprises all components $ = (H a p,V a ), 
the pressure being eliminated by incompressibility. This 
linearisation results in an eigenvalue equation with an op- 
erator £, which acts linearly on the perturbation <fiq(y): 



V'q(y) = £(*(y),<i,dv,dy-)<f><i(y)- 



(4) 



For numerical study, we discretized this equation on a 
staggered mesh. The ID base profile «&(y) was read in 
from the calculation already described. For narrow inter- 
faces, its uniform mesh had too many nodes for use in the 
eigenvalue problem, so we adapted it to put most atten- 
tion near the interface. We then used a NAG routine 
to find the eigenmodes of this discretized problem. 

The results, discussed below, were checked as follows: 
(i) for convergence with respect to mesh structure; (ii) 
that for a homogeneous base state on the underlying con- 
stitutive curve our results match those of Ref. [Rj ; (iii) 
that for a = 0, I = (the local Oldroyd B model), our 
method gives results consistent with Fig. 3 of Ref. [2(| ; 

(iv) that linearisation about a semi-evolved (non-steady) 
banded state using the analytically derived Eqn. 0] gives 
the same results in the limit q x — 0, q z — > as a partic- 
ular direct numerical linearisation performed about the 
same profile in the code that evolves the ID base state; 

(v) for robustness with respect to first evolving the base 
state on either a uniform or adapted grid, using either a 
semi-implicit or explicit algorithm; (vi) that two differ- 
ent methods of eliminating the pressure (using the Oseen 
tensor, and the curl operator) agree. 

For any base profile <&(y) and wavevector q, the num- 
ber of eigenmodes is equal to the number of order pa- 
rameters summed over all mesh points. In this Letter, 
we only consider the eigenvalue w m ax(q) with the largest 
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FIG. 4: Peak of the dispersion relation, i.e. JRtj max at 
dKuWx/dffa = 0. Parameters as for Fig. [3] Limits of the 
banding regime shown by vertical lines. Symbols: data. Dot- 
ted lines: cubic splines, as a guide to the eye. 



real part, 5ftti; max (q). In particular, we ask if this mode is 
stable, 5iw max < 0, or unstable, 5Rw max > 0. All results 
given are for a low solvent viscosity ij = 0.05 -C Gt = 1, 
consistent with experiment. We set a = 0.3, although our 
findings are qualitatively robust to variations in a. This 
leaves the applied shear rate 7 as the tunable parameter. 

The dispersion relation Sw maK fei,9z = 0) for fluctu- 
ations with wavevector confined to the direction of the 
unperturbed flow is shown in Fig. [3] for 7 = 2.0. At 
any q x , 5Rw max increases with decreasing I, and for small 
enough I the dispersion relation is positive over a range 
of wavevectors, rendering the ID profile unstable. For 
small I this applies to shear rates right across the stress 
plateau of Fig. ^ as shown in Fig. 0] Because the I val- 
ues accessed here - I = 0(1 — 10/^m) for a 1mm rheome- 
ter gap - are even larger than those expected physically, 
I = O(100nm), our results suggest that, experimentally, 
the entire stress plateau will be unstable. 

In the limit I — > 0, q x — * 0, the corresponding 
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FIG. 5: Perturbation to flow field si'R-v(y)e iqxX (arrows), 
and contour lines of perturbed normal stress T, xx (x,y) = 
E XI (|/) + S2$lo'xx(y)e t ' ln:X (dotted lines), corresponding to the 
eigenvalue of Fig. |3]with I — 0.01, q x ~ 2.0. Contours down- 
wards: 0.45, 0.60, 0.75, 0.90, 1.05, 1.20, 1.35 (middle value 
shown thicker). Arbitrary scale factors si = 1.5 and S2 = 0.3. 



stress jump AY, XX across the interface (recall Fig|2). This 
triggers a horizontal perturbation to the flow field Qv x in 
these regions, which recirculates, giving an 0(q x ) vertical 



eigenfunction {d y v x , v v = 0, a a p(y)} tends to the spa- velocit y ®°v at ^/2ir = 0.0, 0.5, 1.0. This enhances the 



tial derivative of the base state, d y {d y V x ,V y = 0, S Q /3}, 
representing a simple displacement of the interface in 
the flow-gradient direction, with small corrections in 
the bulk phases to maintain 7 = constant. As q x in- 
creases from zero, this displacement is modulated by a 
wave of wavevector q x yi with an eigenvalue u> max (q x ) — 
+ iq x uj\ + q\^>2 with 0J2 > 0, signifying instability. 
A natural question is whether this instability has the 
same origin as that described by McLeish for the local 
model 16]. It is not obvious, a priori, that this should 
be true because, for the base state at least, the limit 
Z — ^ is singular Q. Indeed, a detailed analysis (work in 
progress) is more complicated in this case, and deferred 
to a longer publication. However, the numerical results of 
Fig. [5] are qualitatively consistent with the mechanism of 
McLeish, as follows. A wavelike interfacial displacement 
with extrema at q x x/2-K = 0.0,0.5,1.0 causes an inter- 
facial tilt near q x x/2n = 0.25,0.75, exposing the normal 



original displacement and so causes instability. Stability 
is restored for higher q x (Fig.[3J), a feature that is absent 
in the local case. 

The eigenvalue 5Rw max (q) over the (q x ,q z ) plane is 
shown in Fig.|SJ Modes with wavevector along the q x axis 
are much more prone to instability than those along the 
q z axis. Nonetheless, for smaller values of I (not shown), 
modes along the line q x = can go unstable as well. 

We note finally an important bound on the validity of 
our calculation. The expansion used to obtain Eqn. 01 
is valid for perturbations that are small at any point in 
space. For example, for the stress components we require 
*C 1. Displacement of the interface by a distance 
A gives cr a /3 = ArfS a( a/dy, which is 0(A/l), because 
the base profile T, a p changes by O(l) over the interfacial 
width 0(1). We are thus restricted to small displace- 
ments, A <ti I. In future work, we will consider A>1. 
We comment briefly on the stability of a sheared in- 
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FIG. 6: Real part of the most unstable eigenvalue, a = 0.3, 
n — 0.05, 7 = 2.0, Reynolds number p/rj — 0.01 (negligible), 
I = 0.01. Contours are -0.45,-0.40... (dotted line), 0.00 
(dashed) and ...0.25,0.30 (solid). 

terface between two phases of a binary fluid in "model 
H" [U. Although this was studied in Ref. 22], that 
work integrated over space to get a simple equation for 
the position of the interface. Such an approach neglects 
changes in the interface's profile, and any fluid flow nor- 
mal to the interface, so is not guaranteed to agree with 
ours. Nonetheless, we found the interface to be stable, as 
in Ref. [13 ■ This supports the idea that normal stresses 
(absent in model H) cause the instability described above. 

To conclude, we have found ID planar shear banded 
flow to be linearly unstable to fluctuations with wavevec- 



tor in the plane of the banding interface, within the non 
local Johnson Segalman model. This applies to shear 
rates right across the stress plateau, suggesting that the 
instability is ubiquitous and therefore that the existing 
theoretical picture of two stable shear bands separated 
by a steady interface needs further thought. Indeed, our 
finding is consistent with accumulating evidence for er- 
ratic fluctuations ^} an d band breakup |2^| in several 
systems. Future work will study the fate of the interface 
in the non linear regime, beyond the validity of this linear 
study. One possibility is that the instability is self limit- 
ing beyond a critical amplitude set by I (e.g., I 1 / 2 ). This 
would be consistent with a narrowly localized but still 
unsteady interface, which might be interpreted as steady 
in experiments that did not have high spatial resolution. 
This might even reconcile early data showing apparently 
steady interfaces with recent work revealing fluctuations. 

By contrast, if the instability were found not to be self 
limiting, and yet ubiquitous in existing banding models 
(work in progress), one would then need a new theoreti- 
cal picture of (reasonably) steady shear bands that could 
still accommodate the required normal stress jump across 
the interface. Other open questions include the status of 
the instability in curved Couette geometry; and the rela- 
tive importance of instabilities at non-zero q (as studied 
here) to those found at zero q in recent models of spatio- 
temporal rheochaos [24|. 
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EPSRC GR/S29560/01 for funding. 



* Electronic address: s.m.fielding@leeds.ac.uk 
[1] M. M. Britton and P. T. Callaghan, Phys. Rev. Lett. 78, 
4930 (1997). 

[2] O. Diat, D. Roux, and F. Nallet, J. Phys. II (France) 3, 
1427 (1993). 

[3] L. Hilliou and D. Vlassopoulos, Ind. Eng. Chem. Res. 41, 
6246 (2002). 

[4] J. S. Raynaud, P. Moucheront, J. C. Baudez, 

F. Bertrand, J. P. Guilbaud, and P. Coussot, J. Rheol. 

46, 709 (2002). 
[5] M. Doi and S. F. Edwards, The Theory of Polymer Dy- 
namics (Clarendon, Oxford, 1989). M. E. Cates, J. Phys. 

Chem. 94, 371 (1990). M. Doi, J. Poly. Sci: Poly. Phys. 

19, 229 (1981). P. D. Olmsted and C.-Y. D. Lu, Phys. 

Rev. E60, 4397 (1999). 
[6] N. A. Spenley, M. E. Cates, and T. C. B. McLeish, Phys. 

Rev. Lett. 71, 939 (1993). P. D. Olmsted and P. M. Gold- 

bart, Phys. Rev. A41, 4578 (1990). 
[7] C.-Y. D. Lu, P. D. Olmsted, and R. C. Ball, Phys. Rev. 

Lett. 84, 642 (2000). 
[8] M. Johnson and D. Segalman, J. Non-Newt. Fl. Mech 2, 

255 (1977). 

[9] P. D. Olmsted, O. Radulescu, and C.-Y. D. Lu, J. Rhe- 
ology 44, 257 (2000). 



[10] A. S. Wunenburger, A. Colin, J. Leng, A. Arneodo, and 
D. Roux, Phys. Rev. Lett. 86, 1374 (2001). R. Bandy- 
opadhyay, G. Basappa, and A. K. Sood, Phys. Rev. Lett. 
84, 2022 (2000). W. M. Holmes, M. R. Lopez-Gonzalez, 
and P. T. Callaghan, Europhys. Lett. 64, 274 (2003). 
Y. T. Hu, P. Boltenhagen, and D. J. Pine, J. Rheol. 42, 
1185 (1998). P. Fischer, E. K. Wheeler, and G. G. Fuller, 
Rheol. Acta 41, 35 (2002). 

[11] S. Manneville, J. B. Salmon, and A. Colin, Eur. Phys. 
J. E 13, 197 (2004). J. B. Salmon, L. Becu, S. Man- 
neville, and A. Colin, Eur. Phys. J. E 10, 209 (2003). 
M. R. Lopez-Gonzalez, W. M. Holmes, P. T. Callaghan, 
and P. J. Photinos, Phys. Rev. Lett. 93, 268302 (2004). 

[12] A. W. El-Kareh and L. G. Leal, J. Non-Newt. Fl. Mech. 
33, 257 (1989). J. L. Goveas and G. H. Fredrickson, Eu- 
rophys. J. B 2, 79 (1998). 

[13] J. Yerushalmi, S. Katz, and R. Shinnar, Chemical Engi- 
neering Science 25, 1891 (1970). 

[14] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical recipes in C (Cambridge University 
Press, Cambridge, 1992). 

[15] Y. Y. Renardy, The. Comp. Fl. Dyn. 7, 463 (1995). 

[16] T. C. B. McLeish, J. Poly. Sci. B-Poly. Phys. 25, 2253 
(1987). 



5 



[17] E. J. Hinch, O. J. Harris, J. M. Rallison, J. Non-Newton. 

Fluid Mech. 43, 311 (1992). 
[18] Numerical Algorithms Group Ltd., Wilkinson House, 

Jordan Hill Road, Oxford, OX2 8DR, UK. 
[19] S. M. Fielding, P. D. Olmsted, Phys. Rev. E 68, 036313 

(2003). 

[20] H. J. Wilson, M. Renardy, and Y. Renardy, J. Non- 
Newton. Fluid Mech. 80, 251 (1999). 

[21] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 
435 (1977). 

[22] A. J. Bray, A. Cavagna, and R. D. M. Travasso, Phys. 

Rev. E 6501, art. no. (2002). 
[23] S Lerouge, PhD Thesis, University of Metz, 2000. 
[24] S. M. Fielding and P. D. Olmsted, Phys. Rev. Lett. 92, 

art. no. (2004). Aradian, A and Cates, M E, Pre-print 

cond-mat/0410509 



